{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeoNeuron/comp-neuro-course/blob/master/tutorials/W0_PythonWorkshop2/W0D2_Tutorial1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "# Python Workshop - LIF Neuron - Part II\n",
    "\n",
    "__Content creators:__ Marco Brigham and the [CCNSS](https://www.ccnss.org/) team\n",
    "\n",
    "__Content Modified:__ Kai Chen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "## Tutorial objectives\n",
    "We learned basic Python and NumPy concepts in the previous tutorial. These new and efficient coding techniques can be applied repeatedly in tutorials from this course, and elsewhere. Please keep them in mind!\n",
    "\n",
    "&nbsp; \n",
    "\n",
    "In this tutorial, we'll introduce spikes in our LIF neuron and evaluate the refractory period's effect in spiking dynamics!\n",
    "\n",
    "&nbsp; \n",
    "\n",
    "Please remember that strengthening your \"Python\" leg will help you advance your \"Neuro\" leg. Alright, let's go for our second walk!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:17.206007Z",
     "iopub.status.busy": "2020-10-15T15:23:17.205143Z",
     "iopub.status.idle": "2020-10-15T15:23:17.489254Z",
     "shell.execute_reply": "2020-10-15T15:23:17.489672Z"
    }
   },
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:17.494081Z",
     "iopub.status.busy": "2020-10-15T15:23:17.493488Z",
     "iopub.status.idle": "2020-10-15T15:23:18.194982Z",
     "shell.execute_reply": "2020-10-15T15:23:18.195580Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Figure settings\n",
    "import ipywidgets as widgets\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:18.208230Z",
     "iopub.status.busy": "2020-10-15T15:23:18.207536Z",
     "iopub.status.idle": "2020-10-15T15:23:18.210097Z",
     "shell.execute_reply": "2020-10-15T15:23:18.209540Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Helper functions\n",
    "\n",
    "t_max = 150e-3   # second\n",
    "dt = 1e-3        # second\n",
    "tau = 20e-3      # second\n",
    "el = -60e-3      # milivolt\n",
    "vr = -70e-3      # milivolt\n",
    "vth = -50e-3     # milivolt\n",
    "r = 100e6        # ohm\n",
    "i_mean = 25e-11  # ampere\n",
    "\n",
    "\n",
    "def plot_all(t_range, v, raster=None, spikes=None, spikes_mean=None):\n",
    "  \"\"\"\n",
    "  Plots Time evolution for\n",
    "  (1) multiple realizations of membrane potential\n",
    "  (2) spikes\n",
    "  (3) mean spike rate (optional)\n",
    "\n",
    "  Args:\n",
    "    t_range (numpy array of floats)\n",
    "        range of time steps for the plots of shape (time steps)\n",
    "\n",
    "    v (numpy array of floats)\n",
    "        membrane potential values of shape (neurons, time steps)\n",
    "\n",
    "    raster (numpy array of floats)\n",
    "        spike raster of shape (neurons, time steps)\n",
    "\n",
    "    spikes (dictionary of lists)\n",
    "        list with spike times indexed by neuron number\n",
    "\n",
    "    spikes_mean (numpy array of floats)\n",
    "        Mean spike rate for spikes as dictionary\n",
    "\n",
    "  Returns:\n",
    "    Nothing.\n",
    "  \"\"\"\n",
    "\n",
    "  v_mean = np.mean(v, axis=0)\n",
    "  fig_w, fig_h = plt.rcParams['figure.figsize']\n",
    "  plt.figure(figsize=(fig_w, 1.5 * fig_h))\n",
    "\n",
    "  ax1 = plt.subplot(3, 1, 1)\n",
    "  for j in range(n):\n",
    "    plt.scatter(t_range, v_n[j], color=\"k\", marker=\".\", alpha=0.01)\n",
    "  plt.plot(t_range, v_mean, 'C1', alpha=0.8, linewidth=3)\n",
    "  plt.xticks([])\n",
    "  plt.ylabel(r'$V_m$ (V)')\n",
    "\n",
    "  if raster is not None:\n",
    "    plt.subplot(3, 1, 2)\n",
    "    spikes_mean = np.mean(raster, axis=0)\n",
    "    plt.imshow(raster, cmap='Greys', origin='lower', aspect='auto')\n",
    "\n",
    "  else:\n",
    "    plt.subplot(3, 1, 2, sharex=ax1)\n",
    "    for j in range(n):\n",
    "      times = np.array(spikes[j])\n",
    "      plt.scatter(times, j * np.ones_like(times), color=\"C0\", marker=\".\", alpha=0.2)\n",
    "\n",
    "  plt.xticks([])\n",
    "  plt.ylabel('neuron')\n",
    "\n",
    "  if spikes_mean is not None:\n",
    "    plt.subplot(3, 1, 3, sharex=ax1)\n",
    "    plt.plot(t_range, spikes_mean)\n",
    "    plt.xlabel('time (s)')\n",
    "    plt.ylabel('rate (Hz)')\n",
    "\n",
    "  plt.tight_layout()\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Sample histogram\n",
    "![histogram statistics](https://github.com/mpbrigham/colaboratory-figures/raw/master/nma/python-for-nma/histogram_statistics.png)\n",
    "\n",
    "Another important statistic is the sample [histogram](https://en.wikipedia.org/wiki/Histogram). For our LIF neuron it provides an approximate representation of the distribution of membrane potential $V_m(t)$ at time $t=t_k\\in[0,t_{max}]$. For $N$ realizations $V\\left(t_k\\right)$ and $J$ bins is given by:\n",
    " \n",
    "\\begin{align*}\n",
    "\\\\\n",
    "N = \\sum_{j=1}^{J} m_j\\\\\n",
    "\\\\\n",
    "\\end{align*}\n",
    "\n",
    "where $m_j$ is a function that counts the number of samples $V\\left(t_k\\right)$ that fall into bin $j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 1\n",
    "Plot an histogram of $J=50$ bins of $N=10000$ realizations of $V(t)$ for $t=t_{max}/10$ and $t=t_{max}$.\n",
    "\n",
    "The function `plt.hist(data, nbins)` plots an histogram of `data` in `nbins` bins. The argument `label` defines a label for `data` and `plt.legend()` adds all labels to the plot.\n",
    "\n",
    "    plt.hist(data, bins, label='my data')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "The parameters `histtype='stepfilled'` and `linewidth=0` may improve histogram appearance (depending on taste). You can read more about [different histtype settings](https://matplotlib.org/gallery/statistics/histogram_histtypes.html).\n",
    "\n",
    "The function `plt.hist` returns the `pdf`, `bins`, and `patches` with the histogram bins, the edges of the bins, and the individual patches used to create the histogram.\n",
    "\n",
    "    pdf, bins, patches = plt.hist(data, bins)\n",
    "\n",
    "We'll make a small correction in the definition of `t_range` to ensure increments of `dt` by using `np.arange` instead of `np.linspace`. (Why is it needed?)\n",
    "\n",
    "    numpy.arange(start, stop, step)\n",
    "\n",
    "**Suggestions**\n",
    "* Plot the 1st histogram after running up to $t=t_{max}/10$\n",
    "* Plot the 2nd histogram at $t=t_{max}$\n",
    "* Add labels to identify the time of each histogram\n",
    "* Plot each histogram with parameters `histtype='stepfilled'` and `linewidth=0`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:18.217109Z",
     "iopub.status.busy": "2020-10-15T15:23:18.216502Z",
     "iopub.status.idle": "2020-10-15T15:23:18.537800Z",
     "shell.execute_reply": "2020-10-15T15:23:18.538162Z"
    },
    "outputId": "bba4acb9-de46-4719-9f7f-0908d89228f0"
   },
   "outputs": [],
   "source": [
    "# set random number generator\n",
    "np.random.seed(2020)\n",
    "\n",
    "# initialize t_range, step_end, n, v_n, syn and nbins\n",
    "t_range = np.arange(0, t_max, dt)\n",
    "step_end = len(t_range)\n",
    "n = 10000\n",
    "v_n = el * np.ones([n, step_end])\n",
    "syn = i_mean * (1 + 0.1 * (t_max / dt)**(0.5) * (2 * np.random.random([n, step_end]) - 1))\n",
    "nbins = 32\n",
    "\n",
    "# loop time steps\n",
    "for step, t in enumerate(t_range):\n",
    "  if step==0:\n",
    "    continue\n",
    "\n",
    "  v_n[:, step] =  v_n[:, step - 1] + (dt / tau) * (el - v_n[:, step - 1] + r * syn[:, step])\n",
    "\n",
    "# initialize the figure\n",
    "plt.figure()\n",
    "plt.ylabel('Frequency')\n",
    "plt.xlabel('$V_m$ (V)')\n",
    "\n",
    "# Complete these lines and uncomment\n",
    "# plt.hist(...)\n",
    "# plt.hist(...)\n",
    "# plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 432
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:18.545042Z",
     "iopub.status.busy": "2020-10-15T15:23:18.544523Z",
     "iopub.status.idle": "2020-10-15T15:23:18.771349Z",
     "shell.execute_reply": "2020-10-15T15:23:18.771852Z"
    },
    "outputId": "7aaa7b56-c4ec-4adc-9753-82b522f68441"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W0_PythonWorkshop/solutions/W0_Tutorial2_Solution_a9562633.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=415 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W0_PythonWorkshop/static/W0_Tutorial2_Solution_a9562633_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Introduce spikes\n",
    "A spike takes place whenever $V(t)$ crosses $V_{th}$. In that case, a spike is recorded and $V(t)$ resets to $V_{reset}$ value. This is summarized in the *reset condition*:\n",
    "$$V(t) = V_{reset}\\quad \\text{ if } V(t)\\geq V_{th}$$\n",
    "\n",
    "For more information about spikes or action potentials see [here](https://en.wikipedia.org/wiki/Action_potential) and [here](https://www.khanacademy.org/test-prep/mcat/organ-systems/neuron-membrane-potentials/a/neuron-action-potentials-the-creation-of-a-brain-signal).\n",
    "\n",
    "\n",
    "![spikes cartoon](https://github.com/mpbrigham/colaboratory-figures/raw/master/nma/python-for-nma/spikes_carton.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 2\n",
    "Insert the reset condition, and collect the spike times of each realization in a dictionary variable `spikes`, with $N=500$.\n",
    "\n",
    "We've used `plt.plot` for plotting lines and also for plotting dots at `(x,y)` coordinates, which is a [scatter plot](https://en.wikipedia.org/wiki/Scatter_plot). From here on, we'll use use `plt.plot` for plotting lines and for scatter plots: `plt.scatter`.\n",
    "```\n",
    "plt.scatter(x, y, color=\"k\", marker=\".\")\n",
    "```\n",
    "\n",
    "A *raster plot* represents spikes from multiple neurons by plotting dots at spike times from neuron `j` at plot height `j`, i.e. \n",
    "```\n",
    "plt.scatter(spike_times, j*np.ones_like(spike_times))\n",
    "```\n",
    "\n",
    "![spikes plot](https://github.com/mpbrigham/colaboratory-figures/raw/master/nma/python-for-nma/spikes_plot.png)\n",
    "\n",
    "In this exercise, we use `plt.subplot` for multiple plots in the same figure. These plots can share the same `x` or `y` axis by specifying the parameter `sharex` or `sharey`. Add `plt.tight_layout()` at the end to automatically adjust subplot parameters to fit the figure area better. Please see the example below for a row of two plots sharing axis `y`.\n",
    "```\n",
    "# initialize the figure\n",
    "plt.figure()\n",
    "\n",
    "# collect axis of 1st figure in ax1\n",
    "ax1 = plt.subplot(1, 2, 1)\n",
    "plt.plot(t_range, my_data_left)\n",
    "plt.ylabel('ylabel')\n",
    "\n",
    "# share axis x with 1st figure\n",
    "plt.subplot(1, 2, 2, sharey=ax1)\n",
    "plt.plot(t_range, my_data_right)\n",
    "\n",
    "# automatically adjust subplot parameters to figure\n",
    "plt.tight_layout()\n",
    "```\n",
    "\n",
    "**Suggestions**\n",
    "* At each time step loop all neurons to:\n",
    "  * Reset $V_n(t)$ to $V_{reset}$ if $V_n(t)\\geq V_{th}$\n",
    "  * Collect spike times in dictionary `spikes` indexed by neuron number\n",
    "  * Add spike count to numpy array `spikes_n`\n",
    "  * In `plt.subplot(3, 1, 2)` loop each neuron to collect and plot spike times\n",
    "  * Plot realizations and sample mean in `plt.subplot(3, 1, 1)`, and mean spiking rate in `plt.subplot(3, 1, 3)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:18.780818Z",
     "iopub.status.busy": "2020-10-15T15:23:18.780275Z",
     "iopub.status.idle": "2020-10-15T15:23:20.432262Z",
     "shell.execute_reply": "2020-10-15T15:23:20.432666Z"
    },
    "outputId": "53243ec6-4d22-4280-bdd2-99e069f80789"
   },
   "outputs": [],
   "source": [
    "# set random number generator\n",
    "np.random.seed(2020)\n",
    "\n",
    "# initialize step_end, t_range, n, v_n and syn\n",
    "t_range = np.arange(0, t_max, dt)\n",
    "step_end = len(t_range)\n",
    "n = 500\n",
    "v_n = el * np.ones([n, step_end])\n",
    "syn = i_mean * (1 + 0.1 * (t_max / dt)**(0.5) * (2 * np.random.random([n, step_end]) - 1))\n",
    "\n",
    "# initialize spikes and spikes_n\n",
    "spikes = {j: [] for j in range(n)}\n",
    "spikes_n = np.zeros([step_end])\n",
    "\n",
    "# loop time steps\n",
    "for step, t in enumerate(t_range):\n",
    "  if step==0:\n",
    "    continue\n",
    "\n",
    "  v_n[:, step] = v_n[:, step - 1] + (dt / tau) * (el - v_n[:, step - 1] + r*syn[:, step])\n",
    "\n",
    "  # collect spike times and reset to resting potential\n",
    "  # Complete these lines and uncomment\n",
    "  # for j in range(n):\n",
    "    # if v_n[j, step] >= vth:\n",
    "      # v_n[j, step] = ...\n",
    "      # spikes[j] += ...\n",
    "      # spikes_n[step] += ...\n",
    "\n",
    "# collect mean Vm and mean spiking rate\n",
    "v_mean = np.mean(v_n, axis=0)\n",
    "spikes_mean =  spikes_n / n\n",
    "\n",
    "# initialize the figure\n",
    "plt.figure()\n",
    "\n",
    "# collect axis of 1st figure in ax1\n",
    "ax1 = plt.subplot(3, 1, 1)\n",
    "for j in range(n):\n",
    "  plt.scatter(t_range, v_n[j], color=\"k\", marker=\".\", alpha=0.01)\n",
    "plt.plot(t_range, v_mean, 'C1', alpha=0.8, linewidth=3)\n",
    "plt.ylabel('$V_m$ (V)')\n",
    "\n",
    "# share xaxis with 1st figure\n",
    "plt.subplot(3, 1, 2, sharex=ax1)\n",
    "# for each neuron j: collect spike times and plot them at height j\n",
    "# Complete these lines and uncomment\n",
    "# for j in range(n):\n",
    "  # times = ...\n",
    "  # plt.scatter(...)\n",
    "\n",
    "plt.ylabel('neuron')\n",
    "\n",
    "# share xaxis with 1st figure\n",
    "plt.subplot(3, 1, 3, sharex=ax1)\n",
    "plt.plot(t_range, spikes_mean)\n",
    "plt.xlabel('time (s)')\n",
    "plt.ylabel('rate (Hz)')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 433
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:20.442677Z",
     "iopub.status.busy": "2020-10-15T15:23:20.442124Z",
     "iopub.status.idle": "2020-10-15T15:23:23.739685Z",
     "shell.execute_reply": "2020-10-15T15:23:23.740032Z"
    },
    "outputId": "1cda77e4-e18f-4f6b-c2d3-0fcdf6923ff2"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W0_PythonWorkshop/solutions/W0_Tutorial2_Solution_a56bd3d0.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=416 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W0_PythonWorkshop/static/W0_Tutorial2_Solution_a56bd3d0_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Boolean indexes\n",
    "Numpy arrays can be indexed by boolean arrays to select a subset of elements (also works with lists of booleans).\n",
    "\n",
    "The boolean array itself can be initiated by a condition, as shown in the example below.\n",
    "```\n",
    "a = np.array([1, 2, 3])\n",
    "b = a>=2\n",
    "print(b)\n",
    "--> [False  True  True]\n",
    "\n",
    "print(a[b])\n",
    "--> [2 3]\n",
    "\n",
    "print(a[a>=2])\n",
    "--> [2 3]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 3\n",
    "We can avoid looping all neurons in each time step by identifying with boolean arrays the indexes of neurons that spiked in the previous step.\n",
    "\n",
    "In the example below, `v_rest` is a boolean numpy array with `True` in each index of `v_n` with value `vr` at time index `step`:\n",
    "\n",
    "    v_rest = (v_n[:,step] == vr)\n",
    "    print(v_n[v_rest,step])\n",
    "    \n",
    "    --> [vr, ..., vr]\n",
    "\n",
    "The function `np.where` returns indexes of boolean arrays with `True` values.\n",
    "\n",
    "You may use the helper function `plot_all` that implements the figure from the previous exercise.\n",
    "\n",
    "**Suggestions**\n",
    "* At each time step:\n",
    "  * Initialize boolean numpy array `spiked` with $V_n(t)\\geq V_{th}$\n",
    "  * Set to `vr` indexes of `v_n` using `spiked`\n",
    "  * Obtain indexes of spiking neurons using `np.where(spiked)[0]`\n",
    "  * Collect spike times in dictionary `spikes` and spike count in numpy array `spikes_n`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "**Helper function(s)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 486
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:23.743616Z",
     "iopub.status.busy": "2020-10-15T15:23:23.743095Z",
     "iopub.status.idle": "2020-10-15T15:23:23.745439Z",
     "shell.execute_reply": "2020-10-15T15:23:23.745852Z"
    },
    "outputId": "63f43a72-4e64-4387-93a4-f5c6599aff8d"
   },
   "outputs": [],
   "source": [
    "help(plot_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 647
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:23.751849Z",
     "iopub.status.busy": "2020-10-15T15:23:23.751347Z",
     "iopub.status.idle": "2020-10-15T15:23:26.759563Z",
     "shell.execute_reply": "2020-10-15T15:23:26.759919Z"
    },
    "outputId": "e57276c1-c64f-4902-cc78-e652073ba151"
   },
   "outputs": [],
   "source": [
    "# set random number generator\n",
    "np.random.seed(2020)\n",
    "\n",
    "# initialize step_end, t_range, n, v_n and syn\n",
    "t_range = np.arange(0, t_max, dt)\n",
    "step_end = len(t_range)\n",
    "n = 500\n",
    "v_n = el * np.ones([n,step_end])\n",
    "syn = i_mean * (1 + 0.1*(t_max/dt)**(0.5)*(2*np.random.random([n,step_end])-1))\n",
    "\n",
    "# initialize spikes and spikes_n\n",
    "spikes = {j: [] for j in range(n)}\n",
    "spikes_n = np.zeros([step_end])\n",
    "\n",
    "# loop time steps\n",
    "for step, t in enumerate(t_range):\n",
    "  if step==0:\n",
    "    continue\n",
    "\n",
    "  v_n[:,step] = v_n[:,step-1] + dt/tau * (el - v_n[:,step-1] + r*syn[:,step])\n",
    "\n",
    "  # collect neurons that spiked and reset to resting potential\n",
    "  spiked = ...\n",
    "  # Complete these lines and uncomment\n",
    "  # v_n[spiked,step] = ...\n",
    "\n",
    "  # collect spike times\n",
    "  # for j in np.where(spiked)[0]:\n",
    "    # spikes[j] += ...\n",
    "    # spikes_n[step] += ...\n",
    "\n",
    "# collect mean spiking rate\n",
    "spikes_mean = spikes_n / n\n",
    "\n",
    "# plot multiple realizations of Vm, spikes and mean spike rate\n",
    "plot_all(t_range, v_n, spikes=spikes, spikes_mean=spikes_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 649
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:26.766808Z",
     "iopub.status.busy": "2020-10-15T15:23:26.766305Z",
     "iopub.status.idle": "2020-10-15T15:23:30.016313Z",
     "shell.execute_reply": "2020-10-15T15:23:30.016841Z"
    },
    "outputId": "d48cc91f-91f1-4168-c2fd-6a0aaa6a6d01"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W0_PythonWorkshop/solutions/W0_Tutorial2_Solution_af1e115b.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=632 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W0_PythonWorkshop/static/W0_Tutorial2_Solution_af1e115b_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 4\n",
    "A *binary raster plot* represents spike times as `1`s in a binary grid initialized with `0`s. We start with a numpy array `raster` of zeros with shape `(neurons, time steps)`, and represent a spike from neuron `5` at time step `20` as `raster(5,20)=1`, for example.\n",
    "\n",
    "The *binary raster plot* is much more efficient than the previous method by plotting the numpy array `raster` as an image:\n",
    "```\n",
    "plt.imshow(raster, cmap='Greys', origin='lower', aspect='auto')\n",
    "```\n",
    "\n",
    "**Suggestions**\n",
    "* At each time step:\n",
    "  * Initialize boolean numpy array `spiked` with $V_n(t)\\geq V_{th}$\n",
    "  * Set to `vr` indexes of `v_n` using `spiked`\n",
    "  * Set to `1` indexes of numpy array `raster` using `spiked`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 646
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:30.026272Z",
     "iopub.status.busy": "2020-10-15T15:23:30.025591Z",
     "iopub.status.idle": "2020-10-15T15:23:31.621042Z",
     "shell.execute_reply": "2020-10-15T15:23:31.621405Z"
    },
    "outputId": "8af8d49e-c133-4f55-daf3-c5eaba28c769"
   },
   "outputs": [],
   "source": [
    "# set random number generator\n",
    "np.random.seed(2020)\n",
    "\n",
    "# initialize step_end, t_range, n, v_n and syn\n",
    "t_range = np.arange(0, t_max, dt)\n",
    "step_end = len(t_range)\n",
    "n = 500\n",
    "v_n = el * np.ones([n,step_end])\n",
    "syn = i_mean * (1 + 0.1*(t_max/dt)**(0.5)*(2*np.random.random([n,step_end])-1))\n",
    "\n",
    "# initialize binary numpy array for raster plot\n",
    "raster = np.zeros([n,step_end])\n",
    "\n",
    "# loop time steps\n",
    "for step, t in enumerate(t_range):\n",
    "  if step==0:\n",
    "    continue\n",
    "\n",
    "  v_n[:,step] = v_n[:,step-1] + dt/tau * (el - v_n[:,step-1] + r*syn[:,step])\n",
    "\n",
    "  # boolean array spiked indexes neurons with v>=vth\n",
    "  spiked = ...\n",
    "  # Complete these lines and uncomment\n",
    "  # v_n[spiked,step] = ...\n",
    "  # raster[spiked,step] = ...\n",
    "\n",
    "# plot multiple realizations of Vm, spikes and mean spike rate\n",
    "plot_all(t_range, v_n, raster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 649
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:31.627503Z",
     "iopub.status.busy": "2020-10-15T15:23:31.626984Z",
     "iopub.status.idle": "2020-10-15T15:23:33.420318Z",
     "shell.execute_reply": "2020-10-15T15:23:33.420679Z"
    },
    "outputId": "41aaa774-0ad5-409c-cb69-85d64cfd9b40"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W0_PythonWorkshop/solutions/W0_Tutorial2_Solution_bc4fc660.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=632 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W0_PythonWorkshop/static/W0_Tutorial2_Solution_bc4fc660_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Refractory period\n",
    "\n",
    "The absolute refractory period is a time interval in the order of a few milliseconds during which synaptic input will not lead to a 2nd spike, no matter how strong. This effect is due to the biophysics of the neuron membrane channels, and you can read more about absolute and relative refractory period [here](https://content.byui.edu/file/a236934c-3c60-4fe9-90aa-d343b3e3a640/1/module5/readings/refractory_periods.html) and [here](https://en.wikipedia.org/wiki/Refractory_period_(physiology)).\n",
    "\n",
    "![refractory period plot](https://github.com/mpbrigham/colaboratory-figures/raw/master/nma/python-for-nma/refractory_period.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 5\n",
    "Investigate the effect of (absolute) refractory period $t_{ref}$ on the evolution of output rate $\\lambda(t)$. Add refractory period $t_{ref}=10$ ms after each spike, during which $V(t)$ is clamped to $V_{reset}$.\n",
    "\n",
    "Experiment increasingly higher values of refractory period $t_{ref}$\n",
    "\n",
    "**Suggestions**\n",
    "* At each time step:\n",
    "  * Manage spike times as in the previous exercise\n",
    "  * Update numpy array `last_spike` with time `t` for spiking neurons\n",
    "  * Initialize boolean numpy array `clamped` using `last_spike`, `t` and `t_ref`\n",
    "  * Reset to `vr` clamped neurons using `clamped`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 646
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:33.427396Z",
     "iopub.status.busy": "2020-10-15T15:23:33.426880Z",
     "iopub.status.idle": "2020-10-15T15:23:35.068315Z",
     "shell.execute_reply": "2020-10-15T15:23:35.068667Z"
    },
    "outputId": "c0defdde-9c4d-4106-a513-dc3311a35367"
   },
   "outputs": [],
   "source": [
    "# set random number generator\n",
    "np.random.seed(2020)\n",
    "\n",
    "# initialize step_end, t_range, n, v_n, syn and raster\n",
    "t_range = np.arange(0, t_max, dt)\n",
    "step_end = len(t_range)\n",
    "n = 500\n",
    "v_n = el * np.ones([n,step_end])\n",
    "syn = i_mean * (1 + 0.1*(t_max/dt)**(0.5)*(2*np.random.random([n,step_end])-1))\n",
    "raster = np.zeros([n,step_end])\n",
    "\n",
    "# initialize t_ref and last_spike\n",
    "t_ref = 0.01\n",
    "last_spike = -t_ref * np.ones([n])\n",
    "\n",
    "# loop time steps\n",
    "for step, t in enumerate(t_range):\n",
    "  if step==0:\n",
    "    continue\n",
    "\n",
    "  v_n[:,step] = v_n[:,step-1] + dt/tau * (el - v_n[:,step-1] + r*syn[:,step])\n",
    "\n",
    "  # boolean array spiked indexes neurons with v>=vth\n",
    "  spiked = (v_n[:,step] >= vth)\n",
    "  v_n[spiked,step] = vr\n",
    "  raster[spiked,step] = 1.\n",
    "\n",
    "  # boolean array clamped indexes refractory neurons\n",
    "  clamped = ...\n",
    "  # Complete these lines and uncomment\n",
    "  # v_n[clamped,step] = ...\n",
    "  # last_spike[spiked] = ...\n",
    "\n",
    "# plot multiple realizations of Vm, spikes and mean spike rate\n",
    "plot_all(t_range, v_n, raster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 649
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:35.156387Z",
     "iopub.status.busy": "2020-10-15T15:23:35.155694Z",
     "iopub.status.idle": "2020-10-15T15:23:36.820380Z",
     "shell.execute_reply": "2020-10-15T15:23:36.820731Z"
    },
    "outputId": "d0fc9a2f-ee97-402e-96b5-aa10d7579c0b"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W0_PythonWorkshop/solutions/W0_Tutorial2_Solution_1ca67e4a.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=632 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W0_PythonWorkshop/static/W0_Tutorial2_Solution_1ca67e4a_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Random refractory period\n",
    "Investigate the effect of random refractory period $t_{ref}$ with \n",
    "$t_{ref} = \\mu + \\sigma\\,\\mathcal{N}$, where $\\mathcal{N}$ is the [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution), $\\mu=0.01$ and $\\sigma=0.007$.\n",
    "\n",
    "Refractory period samples `t_ref` of size `n` is initialized with `np.random.normal`. Remember to clip negative values to `0` with boolean indexes. (Why?)\n",
    "\n",
    "What is the effect of different $\\sigma$ values?\n",
    "\n",
    "**Instructions**\n",
    "* Please execute the cell below\n",
    "* Vary values of $\\sigma$ and observe changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 980,
     "referenced_widgets": [
      "19969f9347c049c9bdd0af96fc65200d",
      "46e98554f2754b648572fe6f273e6468",
      "26960de0547b403988eab42f35a882a4",
      "352d32973d814c59bd85500860bd4b63",
      "c35511ea69724a23a895ea0464f4d300",
      "3c934ca5ded844d19b5bb19d76f692e0",
      "6f44296343b645e59190e2d8272e2e9b",
      "73600740a40b482991b7c2a96e8f75b6",
      "9b7c9e8522f04b8691421a31613e29a2",
      "0e77c726b6a645209d5f6c19043a6220"
     ]
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:36.840797Z",
     "iopub.status.busy": "2020-10-15T15:23:36.840267Z",
     "iopub.status.idle": "2020-10-15T15:23:38.690796Z",
     "shell.execute_reply": "2020-10-15T15:23:38.690346Z"
    },
    "outputId": "da29db82-7a89-4eeb-cacf-5547abbbf32c"
   },
   "outputs": [],
   "source": [
    "# @title Interactive demo: Parameter exploration\n",
    "\n",
    "def random_ref_period(mu, sigma):\n",
    "  # set random number generator\n",
    "  np.random.seed(2020)\n",
    "\n",
    "  # initialize step_end, t_range, n, v_n, syn and raster\n",
    "  t_range = np.arange(0, t_max, dt)\n",
    "  step_end = len(t_range)\n",
    "  n = 500\n",
    "  v_n = el * np.ones([n,step_end])\n",
    "  syn = i_mean * (1 + 0.1*(t_max/dt)**(0.5)*(2*np.random.random([n,step_end])-1))\n",
    "  raster = np.zeros([n,step_end])\n",
    "\n",
    "  # initialize t_ref and last_spike\n",
    "  t_ref = mu + sigma*np.random.normal(size=n)\n",
    "  t_ref[t_ref<0] = 0\n",
    "  last_spike = -t_ref * np.ones([n])\n",
    "\n",
    "  # loop time steps\n",
    "  for step, t in enumerate(t_range):\n",
    "    if step==0:\n",
    "      continue\n",
    "\n",
    "    v_n[:,step] = v_n[:,step-1] + dt/tau * (el - v_n[:,step-1] + r*syn[:,step])\n",
    "\n",
    "    # boolean array spiked indexes neurons with v>=vth\n",
    "    spiked = (v_n[:,step] >= vth)\n",
    "    v_n[spiked,step] = vr\n",
    "    raster[spiked,step] = 1.\n",
    "\n",
    "    # boolean array clamped indexes refractory neurons\n",
    "    clamped = (last_spike + t_ref > t)\n",
    "    v_n[clamped,step] = vr\n",
    "    last_spike[spiked] = t\n",
    "\n",
    "  # plot multiple realizations of Vm, spikes and mean spike rate\n",
    "  plot_all(t_range, v_n, raster)\n",
    "\n",
    "  # plot histogram of t_ref\n",
    "  plt.figure(figsize=(8,4))\n",
    "  plt.hist(t_ref, bins=32, histtype='stepfilled', linewidth=0, color='C1')\n",
    "  plt.xlabel(r'$t_{ref}$ (s)')\n",
    "  plt.ylabel('count')\n",
    "  plt.tight_layout()\n",
    "\n",
    "_ = widgets.interact(random_ref_period, mu = (0.01, 0.05, 0.01), \\\n",
    "                              sigma = (0.001, 0.02, 0.001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Using functions\n",
    "Running key parts of your code inside functions improves your coding narrative by making it clearer and more flexible to future changes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 6\n",
    "Re-organize part of the code from the previous exercise with functions.\n",
    "\n",
    "**Suggestions**\n",
    "* Complete the function `spike_clamp()` to update $V(t)$ and deal with spiking and refractoriness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 646
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:38.708273Z",
     "iopub.status.busy": "2020-10-15T15:23:38.707721Z",
     "iopub.status.idle": "2020-10-15T15:23:40.299390Z",
     "shell.execute_reply": "2020-10-15T15:23:40.299737Z"
    },
    "outputId": "131d71fe-efe7-47e3-dd1e-4246e7103090"
   },
   "outputs": [],
   "source": [
    "def ode_step(v, i, dt):\n",
    "  \"\"\"\n",
    "  Evolves membrane potential by one step of discrete time integration\n",
    "\n",
    "  Args:\n",
    "    v (numpy array of floats)\n",
    "      membrane potential at previous time step of shape (neurons)\n",
    "\n",
    "    v (numpy array of floats)\n",
    "      synaptic input at current time step of shape (neurons)\n",
    "\n",
    "    dt (float)\n",
    "      time step increment\n",
    "\n",
    "  Returns:\n",
    "    v (numpy array of floats)\n",
    "      membrane potential at current time step of shape (neurons)\n",
    "  \"\"\"\n",
    "  v = v + dt/tau * (el - v + r*i)\n",
    "\n",
    "  return v\n",
    "\n",
    "def spike_clamp(v, delta_spike):\n",
    "  \"\"\"\n",
    "  Resets membrane potential of neurons if v>= vth\n",
    "  and clamps to vr if interval of time since last spike < t_ref\n",
    "\n",
    "  Args:\n",
    "    v (numpy array of floats)\n",
    "      membrane potential of shape (neurons)\n",
    "\n",
    "    delta_spike (numpy array of floats)\n",
    "      interval of time since last spike of shape (neurons)\n",
    "\n",
    "  Returns:\n",
    "    v (numpy array of floats)\n",
    "      membrane potential of shape (neurons)\n",
    "    spiked (numpy array of floats)\n",
    "      boolean array of neurons that spiked  of shape (neurons)\n",
    "  \"\"\"\n",
    "\n",
    "  ####################################################\n",
    "  ## TODO for students: compute the mean squared error\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError(\"Student exercise: complete the function\")\n",
    "  ####################################################\n",
    "  # boolean array spiked indexes neurons with v>=vth\n",
    "  spiked = ...\n",
    "  v[spiked] = ...\n",
    "\n",
    "  # boolean array clamped indexes refractory neurons\n",
    "  clamped = ...\n",
    "  v[clamped] = ...\n",
    "\n",
    "  return v, spiked\n",
    "\n",
    "# set random number generator\n",
    "np.random.seed(2020)\n",
    "\n",
    "# initialize step_end, t_range, n, v_n, syn and raster\n",
    "t_range = np.arange(0, t_max, dt)\n",
    "step_end = len(t_range)\n",
    "n = 500\n",
    "v_n = el * np.ones([n,step_end])\n",
    "syn = i_mean * (1 + 0.1*(t_max/dt)**(0.5)*(2*np.random.random([n,step_end])-1))\n",
    "raster = np.zeros([n,step_end])\n",
    "\n",
    "# initialize t_ref and last_spike\n",
    "mu = 0.01\n",
    "sigma = 0.007\n",
    "t_ref = mu + sigma*np.random.normal(size=n)\n",
    "t_ref[t_ref<0] = 0\n",
    "\n",
    "# initialize t_ref and last_spike\n",
    "last_spike = -t_ref * np.ones([n])\n",
    "\n",
    "# loop time steps\n",
    "for step, t in enumerate(t_range):\n",
    "  if step==0:\n",
    "    continue\n",
    "\n",
    "  v_n[:,step] = ode_step(v_n[:,step-1], syn[:,step], dt)\n",
    "  # Complete the`spike_clamp` function above and uncomment\n",
    "  # v_n[:,step], spiked = spike_clamp(v_n[:,step], t - last_spike)\n",
    "\n",
    "  # update raster and last_spike\n",
    "  # raster[spiked,step] = 1.\n",
    "  # last_spike[spiked] = t\n",
    "\n",
    "# plot multiple realizations of Vm, spikes and mean spike rate\n",
    "plot_all(t_range, v_n, raster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 649
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:40.390346Z",
     "iopub.status.busy": "2020-10-15T15:23:40.389636Z",
     "iopub.status.idle": "2020-10-15T15:23:42.066478Z",
     "shell.execute_reply": "2020-10-15T15:23:42.066995Z"
    },
    "outputId": "301d098e-466c-42da-d4b9-687d4f4ed419"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W0_PythonWorkshop/solutions/W0_Tutorial2_Solution_80589e38.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=632 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W0_PythonWorkshop/static/W0_Tutorial2_Solution_80589e38_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Using classes\n",
    "Using classes helps with code reuse and reliability. Well-designed classes are like black boxes in that they receive inputs and provide expected outputs. The details of how the class processes inputs and produces outputs are unimportant.\n",
    "\n",
    "See additional details here: [A Beginner's Python Tutorial/Classes](https://en.wikibooks.org/wiki/A_Beginner%27s_Python_Tutorial/Classes)\n",
    "\n",
    "*Attributes* are variables internal to the class, and *methods* are functions internal to the class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 7\n",
    "In this exercise we'll practice with Python classes by implementing `LIFNeurons`, a class that evolves and keeps state of multiple realizations of LIF neurons.\n",
    "\n",
    "Several attributes are used to keep state of our neurons:\n",
    "```\n",
    "self.v             current membrane potential\n",
    "self.spiked        neurons that spiked\n",
    "self.last_spike    last spike time of each neuron\n",
    "self.t             running time of the simulation\n",
    "self.steps         simulation step\n",
    "```\n",
    "There is a single method:\n",
    "```\n",
    "self.ode_step()    performs single step discrete time integration\n",
    "                   for provided synaptic current and dt\n",
    "```\n",
    "\n",
    "**Suggestions**\n",
    "* Complete the spike and clamp part of method `self.ode_step` (should be similar to function `spike_and_clamp` seen before)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 664
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:42.078299Z",
     "iopub.status.busy": "2020-10-15T15:23:42.077771Z",
     "iopub.status.idle": "2020-10-15T15:23:43.712423Z",
     "shell.execute_reply": "2020-10-15T15:23:43.712775Z"
    },
    "outputId": "b6324e62-a195-4e1d-ef8c-eac70ddae50f"
   },
   "outputs": [],
   "source": [
    "# simulation class\n",
    "class LIFNeurons:\n",
    "  \"\"\"\n",
    "  Keeps track of membrane potential for multiple realizations of LIF neuron,\n",
    "  and performs single step discrete time integration.\n",
    "  \"\"\"\n",
    "  def __init__(self, n, t_ref_mu=0.01, t_ref_sigma=0.002,\n",
    "               tau=20e-3, el=-60e-3, vr=-70e-3, vth=-50e-3, r=100e6):\n",
    "\n",
    "    # neuron count\n",
    "    self.n = n\n",
    "\n",
    "    # neuron parameters\n",
    "    self.tau = tau        # second\n",
    "    self.el = el          # milivolt\n",
    "    self.vr = vr          # milivolt\n",
    "    self.vth = vth        # milivolt\n",
    "    self.r = r            # ohm\n",
    "\n",
    "    # initializes refractory period distribution\n",
    "    self.t_ref_mu = t_ref_mu\n",
    "    self.t_ref_sigma = t_ref_sigma\n",
    "    self.t_ref = self.t_ref_mu + self.t_ref_sigma * np.random.normal(size=self.n)\n",
    "    self.t_ref[self.t_ref<0] = 0\n",
    "\n",
    "    # state variables\n",
    "    self.v = self.el * np.ones(self.n)\n",
    "    self.spiked = self.v >= self.vth\n",
    "    self.last_spike = -self.t_ref * np.ones([self.n])\n",
    "    self.t = 0.\n",
    "    self.steps = 0\n",
    "\n",
    "\n",
    "  def ode_step(self, dt, i):\n",
    "\n",
    "    # update running time and steps\n",
    "    self.t += dt\n",
    "    self.steps += 1\n",
    "\n",
    "    # one step of discrete time integration of dt\n",
    "    self.v = self.v + dt / self.tau * (self.el - self.v + self.r * i)\n",
    "    ####################################################\n",
    "    ## TODO for students: complete the `ode_step` method\n",
    "    # Fill out function and remove\n",
    "    raise NotImplementedError(\"Student exercise: complete the method\")\n",
    "    ####################################################\n",
    "    # spike and clamp\n",
    "    self.spiked = ...\n",
    "    self.v[self.spiked] = ...\n",
    "    self.last_spike[self.spiked] = ...\n",
    "    clamped = ...\n",
    "    self.v[clamped] = ...\n",
    "\n",
    "    self.last_spike[self.spiked] = ...\n",
    "\n",
    "\n",
    "# set random number generator\n",
    "np.random.seed(2020)\n",
    "\n",
    "# initialize step_end, t_range, n, v_n, syn and raster\n",
    "t_range = np.arange(0, t_max, dt)\n",
    "step_end = len(t_range)\n",
    "n = 500\n",
    "v_n = el * np.ones([n,step_end])\n",
    "syn = i_mean * (1 + 0.1*(t_max/dt)**(0.5)*(2*np.random.random([n,step_end])-1))\n",
    "raster = np.zeros([n,step_end])\n",
    "\n",
    "# initialize neurons\n",
    "neurons = LIFNeurons(n)\n",
    "\n",
    "# loop time steps\n",
    "for step, t in enumerate(t_range):\n",
    "  # Complete the method and uncomment\n",
    "  # neurons.ode_step(dt, syn[:,step])\n",
    "\n",
    "  # log v_n and spike history\n",
    "  v_n[:,step] = neurons.v\n",
    "  raster[neurons.spiked, step] = 1.\n",
    "\n",
    "# report running time and steps\n",
    "print(f'Ran for {neurons.t:.3}s in {neurons.steps} steps.')\n",
    "\n",
    "# plot multiple realizations of Vm, spikes and mean spike rate\n",
    "plot_all(t_range, v_n, raster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 666
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2020-10-15T15:23:43.723880Z",
     "iopub.status.busy": "2020-10-15T15:23:43.723364Z",
     "iopub.status.idle": "2020-10-15T15:23:45.419907Z",
     "shell.execute_reply": "2020-10-15T15:23:45.420259Z"
    },
    "outputId": "2c63a728-625a-4249-dd2a-62b49a00b428"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W0_PythonWorkshop/solutions/W0_Tutorial2_Solution_86dc40ee.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=632 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W0_PythonWorkshop/static/W0_Tutorial2_Solution_86dc40ee_1.png>\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W0_Tutorial2",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
